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Mesh motion is the process by which a computational domain is updated in time to 
reflect physical changes in the material the domain represents. Such a technique is needed 
in the study of the thermal response of ablative materials, which erode when strong heating 
is applied to the boundary. Traditionally, the thermal solver is coupled with a linear elastic 
or biharmonic system whose sole purpose is to update mesh node locations in response 
to altering boundary heating. Simple mesh motion algorithms rely on boundary surface 
normals. In such schemes, evolution in time will eventually cause the mesh to intersect and 
“tangle” with itself, causing failure. Furthermore, such schemes are greatly limited in the 
problems geometries on which they will be successful. This paper presents a comprehen- 
sive and sophisticated scheme that tailors the directions of motion based on context. By 
choosing directions for each node smartly, the inevitable tangle can be completely avoided 
and mesh motion on complex geometries can be modeled accurately. 


I. Introduction 


Problems involving mesh motion—which are not to be confused with moving mesh methods, a class of 
adaptive mesh redistribution techniques designed to increase solution accuracy at low cost—arise naturally 
in the study of the thermal response of ablative materials. Ablation is the process by which a surface erodes 
by conversion of mass to another state of matter, most often from solid to gas (sublimination). While mass 
is conserved when the entire system is considered, from the perspective of the solid, the mass loss is absolute: 
no residue or accumulation is left, such as when a candle melts. Because the shape of the material changes 
in time, so does the computational domain. This is the goal of mesh motion: independent of the thermal 
solver, update the computational mesh to model its evolution in time. 


A detailed discussion of the physics and modeling of ablation is beyond the scope of this paper. Instead, we 
assume that the total mass loss Am is known. This quantity is related to the mass loss rate due to ablation 
m in the following manner: 


ts 
Am = i: m dS dt (1) 
to an 


By solving a surface energy balance equation, the rate m can easily be obtained at surface quadrature points. 
Assuming a constant material density p, 7 can be converted to the surface recession rate 5: 


s=— 
p 


[5] = m/s [rr] = ke/m-s? 


Because the thermal response problem from which Am is solved at discrete time steps, the time integral in 
(1) must be approximated. This paper assumes a fully explicit approach: the surface recession, determined 
at time level n, is taken constant until time level n + 1. Under this assumption, the total change in surface 
position over a period of At = tn1i1 — ty is 
As = sAt 
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With this quantity known, we can determine the overall mesh motion by assuming the solid responds to 
the strain of deformation by As linearly. The dynamics of displacement u of a position within the material 
(domain) 9 are governed by the linear elasticity equation 


Lu=f xEQ 
(2) 


u = (sAt)d(x) x€ dQ 
where u(-) € R@ and the operator £ is given by 
Lust -VX(V-u)—(V-pV)u— V- (Vu)? 


In general, the homogenous case f = 0 is the only one of practical interest. The vector d indicates the 
direction of the surface motion. The variational problem corresponding to (2) is 


a(u,v) = (£,4) :2(0:R4) Vv € H'(0;R¢) 


where the bilinear form is 


a(u,v) = | (V -u)(V-v) dx +n f (Vu, Vv) paxa OX +n f (Vu, VV) paxa AX 
a Q 


: (3) 
bay -u)(v, V) pa dS — ef (Vu )d,V) pa dS — of ((Vu)d,v) 94 ds 


where we use the standard inner products 


d d 
(A, By igs = trace(A7 B) = ye S° Aj By; 


i=1 j=1 
d 
(£,8) p2@,R9) = : (£,8) pa dx = SS | fe9a00 dx 


Because the sole boundary condition in the problem is of Dirichlet type and will be enforced as an essential 
boundary condition, we require v € H4(Q;R7) so that the surface integrals of (3) all vanish. 


We distinguish three disjoint subsets of OQ: 


OQR: Receding portion Specified motion of As along d 
OQs: Sliding portion Motion in the normal direction forbidden 
OQ: Stationary (fixed) portion No motion in any direction permitted 


In a simple approach, the direction of motion d is taken to be the surface normal v. While simple to 
implement, this has two unfortunate and undesirable consequences: 


1. As we shall see in the coming section, it precludes sliding motion along any axis other than a coordinate 
one in a finite-element based formulation of (2). 


2. On curved surfaces, the normal vectors from element to element are not parallel. This dooms the mesh 
to an inevitable collision with itself. An example is depicted in Figure 1 below. 


This report presents a scheme for two- and three-dimensional meshes that simulatenously removes the first 
restriction while not suffering from the second drawback. We dub this method “tangle-free mesh motion.” 


In one-dimension, there is only one possible direction of motion: “normal” (left or right), always directed 
into the material. Tangling is impossible and sliding conditions do not exist, since prohibiting motion in the 
normal is equivalent to quashing all motion. As a result, mesh motion in one-dimension is fully resolved. 
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Fig. 1: Collision Course. The starting geometry is depicted on the left in (a); here, OQR is the curved 

part and the two flat sides (aligned with the x and y axes) comprise 00g. When aerodynamic heating is 

applied to OQ R of this ablating isoq, after sufficient simulation time, the result is (b). The material has 
receded and the mesh has become tangled. The thermal solve will fail shortly after this occurs. 


II. Two-Dimensional Meshes 


A. General Framework 


Suppose 2 has been discretized with a finite element mesh YT and let 7 € T be an element in this mesh such 
that |T 1 OQ|pa-1 > 0; that is, the intersection with the boundary is nontrivial (more than a single point in 
two-dimensions). Let {x1,--- ,Xm}C OTNOQ and assume that {xXn41,-.-,xv}NOQ =. On T, the finite 
element solution to (2) is a given by a representation in terms of local shape functions: 


N 
ul, = >> ante (x)er + Brbe(x)er (4) 
k=1 


As defined in finite element software packages (e.g., libMesh, deal.II), a shape function is the Lagrange 
interpolant at the nodes of a single element 7*. It has the property w(x,;) = 6;; whenever x, is a node of J*. 
The finite element basis functions are assembled from piecewise definitions involving these shape functions. 


While representing vector solutions with scalar basis functions by decomposing into components as in (4) 
is advantageous in software, it is inconvenient for this discussion. We thus consider a reindexing scheme by 
creating vector-valued shape functions 


7 Wr (x) e1 1<k<N 
PEO) =) ay. w(xlep N+ISk<2N 


so that u is more concisely written as 


2N 
ul. = S 7 ange (x) (5) 
k=1 


It is this separation into Cartesian components that prohibits sliding along any axis other than a coordinate 
one. The boundary condition is a Dirichlet one of movement specified along d. Suppose we enforce this 
condition by setting the appropriate degrees of freedom in the above representation. In R?, if d 4 ae; for 
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some @ and i, then for each k such that x, € OQ, we must set the boundary conditions 


2 5(x,)At =. $(x,_n )At 
RS ae) RNS =a) 


whereby d“ represents the i*” component* of d, all of which are necessarily nonzero by assumption. On a 
sliding set, 6 = 0. All available degrees of freedom go to representing d in the standard basis so that the 
above quashes all motion when d is not aligned with a coordinate axis. 


The solution, therefore, is to change the basis from the standard fer, eo} to one more appropriate for local 
conditions. Define, instead, the following vector basis functions: 


We(X) De l<k<m 
3 We (x) e1 mt1l<k<N 
Pi (x) = 
We—Nn(X) Gk—-Nn N+1<k<m+N 
Wr—n(X) €2 m+N+1<k<2N 





For each k, we require R? = span {px, qx} so that all solution vectors are representable in {@,})_,. For the 
moment, suppose that each px and qx is known a priori; defining these directions is a matter of circumstance 
and the focus of a future section. The new shape functions can be written in terms of the old: 











Pr p, (x) + pp, w(x) l<k<m 
2 Px (X) m+t1l<k<N 
Pu(X) = 4 aa) (2) N+1<k<m+N 6) 
INP K—N(X) + Uw PK(X) es ear 
p, (x) m+N+1<k<2N 
We may encode the transformation as matrix multiplication 
2N 
Pr(x) = >> Muyy(x) (7) 
l=1 
where M,, are the entries of the matrix 
pl 0 0 LQ p? 0 0: 0 
0 ee 0 0 ae 0.0 0. |) 
: | 
0 G00 1] 0 0 0 0 || 
q 9.8 0) G20 0 
0 gio 0| 9 g@®io 0 e 
0 0°: 0 0 : ° 1 0 s 
Co ee ae ; 
0 i= OF oe 0" |I0 0 0 1 || 
UN." ————" 
N-m N-m 
———__——_’ e———__—_———~ 
N N 


The element stiffness matrix K is the Gram matrix for the associated bilinear form a(-,-) when the basis 
functions are restricted only to a single element. It has entries 


Ki = a(%;,9;) 





*We follow standard convention: boldface for vectors quantities, normal weight for scalar ones. 
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Let K be the element stiffness matrix for the new basis functions ~,-. Then using (7), we have 


Ki = a(;, P;) 


2N 2N 
=a (>: Miver(x), > Mju ey «) 
l=1 ah 
2N 2N 


= SO Se Myva(ey(x), Gy (x) Mit 


l=1 V=1 
2N /2N 
= S- (>: Mivku Mii 
i=1 \V=1 
which we recognize, though matrix multiplication, as the congruence relation 
K = M'KM (8) 

The element contribution to righthand-side vector, which we denote as the vector E which has entries 

Ex a (f, Pr) 12(7) 
can be transformed in a similar manner. Let E have entries 

Ey, = (£, Pe) 12(7) 


Then using (7), we have 
2N 


Ey — > Mui(£, 1) 12(7) 
i=1 
or in terms of matrix-vector multiplication 5 
E=M’'E (9) 


A bit of forward thinking reveals that that one-dimensional subspaces are sufficient to describe the motion 
of a node. Once we’ve determined the new location of a receding node, we can compute a single vector of 
motion by subtracting the new position from the old. For sliding sets, motion is forbidden in the normal 
direction but allowed in the orthogonal complement of the normal—this is again a one-dimensional subspace. 
Consequently, nonzero motion will occur in the direction of at most one of px and qx. We can therefore 
simplify matters considerably by choosing that q places M into a favorable form. By taking 


ax = [p?,-p?]" (10) 


M becomes symmetric: 


po? g|o --- 0 p 0 0 LO 











If |p,| = 1, then M becomes an orthogonal matrix, making inversion trivial. 
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B. Applying the Corrections 


As long as M is invertible, we can obtain the stiffness for the old basis vectors from that of the new ones. This 
is especially convenient in software because the component-decomposition of y, greatly simplifies storage 
and bookeeping requirements. For a true representation involving the @;, we must store thousands of vector 
components and track what maps where. It is therefore more advantageous to use the ¢ as an intermediate 
local step and ultimately return to a representation in y,. 


As long as M is invertible, we can do exactly what we propose above by reversing (8) and (9): 
K=M KM! =MKM 
E=M-TE=ME 

with the second equalities following from the orthorgonality of M. 


The symmetric M was constructed under the assertation that only a single direction would describe the 
motion, but this does not mean that we will always enforce a single direction of motion. When only one 
direction is enforced, the dynamics of the linear elastic system determine the motion in the reciprocal basis 
direction. This is precisely what we desire for sliding side sets but not receding ones, on which we wish to 
specify exactly the magnitude and direction of motion. 


Analagous to (5), we can write the finite element solution in terms of the new basis as 
2N 
ul = > Gn Gi(%) 
k=1 


Choose px = d(x,), the direction of the boundary condition at x,. Then we can set the boundary condition 
exactly by equating this condition with the corresponding degree of freedom: 


Qk = 3(Xx) 


where § is the (perhaps nonphysical) prescribed value of recesssion, and when appropriate (i-e., motion in 
the direction of q, should be quashed) 


Apk—n = 0 


Ordinarily, when degrees of freedom are set in this manner, the corresponding variables are removed from 
the system and references to those quantities replaced by the predetermined value. Unfortunately, this is 
not possible here because the above sets degrees of freedom in #;,, which are local shape functions, not finite 
element basis functions. Accordingly, removal of system variables would require a complex and cumbersome 
bookeeping step that would negate much of the perceived value of reducing the system size. 


Instead, we exploit the nature of finite precision arithmetic: if we have 


e lat+tb=e le 


then a + c as ¢ + 0+ if min {|¢|,|2|} < e~!. On a computer, by taking e~! ~ 101°, for instance, it is 
possible to obtain a = c to nearly machine precision. This technique is known as a penalty method. 








To enforce the condition via the penalty method, we modify K and E by adding perturbations: 
Ki+K+P; 
EVE+ B; 

The subscripts on P and B match the number of degrees of freedom being set (single or double enforcement). 


With these modifications, the contributions from the element 7 to the global stiffness matrix A and right- 
hand vector F become, per (8), 


A|-=KOM(K+P))M = Al|,=K+ 
F|,=EGM(E+B) = F\,=E-+ 


MPM 
(12) 
MB 


» 
te *» 
[I> 





A 


= 
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1. Single Enforcement 


When only one degree of freedom is set, the perturbations take the form 


5(x1) 
1/0 hs me 
P, 4-7} B, 4.7! a (13) 
E 
oho: : 
0 


in which I denotes the identity matrix. If we partition M as 


M= Miu M2 
Mo M22 


with the blocks defined along the thick dividing lines of (11), we can exploit the block structure of P and 
perform the required multiplication easily: 


(14) 


MPM = M?, MiiMie2 
MoaiMi MaiMie2 


We thus obtain a nicely simplifed result for the corrections to the stiffness matrix and righthand-side vector: 





a") Q | ao 0 oe 
[pb] pe? im) a? 
0 0 0 

& 

| 

Road () 0 0 0 ge _i| ° a 
TP 0| eT o| feu 
Dp? [| 5 (mdi 

0 0 0 : 

| 

0 0 0 0 0 = 





2. Double Enforcement 


When both degrees of freedom are set, we have 


pagt| & 2-0-2 Beet. oe (15) 
001 0 S(PT,---;Pin) 
0000 0 


The block sizes in Pz mirror those of M when partitioned along the dotted and thin black lines in (11): 


Mi O Mi3 O 
Mie 0 I 0 0 
M3; 0 M33 O 
0 0 0 I 
7 of 30 


American Institute of Aeronautics and Astronautics 


The vector py is q, as defined by (10). Within Bo, the block 8(di,...,dm) € R” is defined as 
§(x1;d1) 
3(Xm3 dm) 
x;) is in the direction d;. Performing the 


The inclusion of a direction d; indicates that the recession value §( 
multiplication to form Ky, we have 


Mj, + M13M31 0 MiiMi3 + Mi3M33 0 


K 0 0 0 0 
2 — 

M31Mi; + M33M3; 0 M31Mj3 + M3; 0 

0 0 0 0 


Each block Mj, is a diagonal matrix. Letting [A] ij denote the (7,7) entry of the matrix A, we can compute 


the entries in the above representation: 


«=P + PPP =3 


ot! Dy — ppp) = 0 


[MiiMi3 + Mi3Ms33| 
ls (1) (1), (2) _ =H 
3] 


[M31Mii + M33M31 =p? dp — Pi Pj 


«= PP] + bP =1 


Then, we obtain an amazingly simple result for K» and slightly complicated one for E: 


[Ms:Mi3 + M3 


(pi,8 (x1; Pi, Pr )) 








1 0 
; : Pm,S Sta Pm; Pn) 
1 0 E = 0 
eincensdnadoencrnearagns 0 ondriuarnds| sre orrneree hee | g 
K = 1 : 0 : 0 gE 1 0 < 
2 7” : ; 2 > 
E10 : 1 : E (pr. 8 (5 Pa, Pr )) 
= 
tee cuascellcbeae. De tee, | (Pin 8(Xm; Pm; Pin) 
- 5 ‘ : g 
0 Bp | 
m m 0 
ee eS “ 
N-m N-m 


We reuse the previous notation (16) for compactness; for clarity, the entries of E, are of the form 
(2) 
(1) 


(p, 8(x*; p, p*)) = 8(x*; p)p™ + 8(x*;p*)p 


(p*,8(x*;p, p+)) = 8(x*; p)p®) — 8(x*; pt)p 








C. Setting the Direction and Magnitude of Motion 


The previous section supposes that the recession values and directions of motion are known a priori. This 
section describes how to determine them. First, we need to define a simple but critical concept: 
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Definition (Side set). 


A side set is a subset of OQ. Each of 0OR, OQg, and OOF is the union of side sets. 





In implementation, the user will specify the side sets so that OQr, OQOg, and OOF are fully determined. For 
any two side sets S and S’, it is not the case that SMS’ = ; instead, we should have the property that 
|S S'|pa-1 = 0 so that non-empty intersections can be no more than singletons in 2D and lines in 3D. 


How the directions of motion are constructed depend on the problem geometry and the location of the node. 
Side sets allow for quick and easy classification that eliminates the need for complex geometric detection 
algorithms. With their help, we reduce the analysis to two node types in two-dimensions: 





Definition (2D Node Classification). 





In two-dimensional meshes, a node is a corner node if it is contained in two side sets. It is an 
interior node if it is contained in only one side set. 





It is important to note the distinction between our definition of a corner and the geometric notion of a 
corner. Automatic detection based on the latter is precarious, as it is easy for an algorithm to misclassify 
coarse portions of the mesh, which can occur even in well-constructed meshes after extensive recession. 


1. Corner Nodes 


In two dimensions, a node is contained in exactly two boundary faces (edge). It is possible for a node to 
be contained in an arbitrary number of elements but our requirement that intersection with the boundary 
be non-trivial excludes from consideration all but one or two. It is beneficial to realize that the elements 
themselves are not important—it is the edges that give us the important information. We therefore simplify 
our analysis by disregarding the elements once we’ve identified the boundary edges. 


The physical recession values § are known only at quadrature points along the edge, yet the results of the 
previous section require a value at the node itself. A potential solution is apparent: 

1. Choose the value of s at the closest quadrature point for each edge. 

2. Interpolate between these two values to construct an estimate for $ at the node. 


While seemingly sound, this method has a major deficiency: with only two values, we must use linear 
interpolation in step 2. Besides physical incorrectness, this effort is unable to consider curvature or other 
properties of the mesh and is likely to lead to poor mesh quality and skewed geometry. Instead, the solution 
is, in essence, to move the entire edge at once: 


1. For each face, compute the “new quadrature points” by moving s along v, both of which are known 
exactly, and then compute the best fit line through these new points to determine the new edge. 


2. Compute the intersection of the new edges and make that the new location of the node. 
3. The motion § and direction vector p are obtained by subtracting the old location from the new. 


This idea is realized in Algorithm 1 below. Further explanation of the steps follows the presentation. 





Algorithm 1 (Corner Motion). 





. Let € be the edge containing the node x* from one side set and €’ be the other. 


. Let x,,...,Xq be the quadrature points of €. Let v, and 8, be the normal and recession 
values, respectively, specified at the quadrature point x,. Let x4,... Xr v',, and $4, be the 
similarly named and corresponding quantities from €’. 
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continued). 


3. Compute the centroids of the sets {x, + (s,At)vn}i_, and {xj + (3, Atv), bi: 
q 
cs S- [xk + (8, At)v¢] 
k=1 
. Form the 2 x q matrix A and 2 x q’ matrix A’ 


A 2 | oat (& At —e | +» | xq 








A’ | xh + (At, —e! | + | x4 


. Compute the singular value decompositions of A and A’: 
A=UZV A’=U'D’V' 


where U € R22, 5 € R2*7, V € R24, U! € R2X2, 5’ € R2*V, and VW’ € RY*%. Let 
ni = ug and ny = us, the third columns of U and U’, respectively. 


a, ax 7 =c’-c 
Ny —ns t* 
and set x*,,, =ce+s*[—- ni), nr]. 


. Solve the linear system 


—x*| and p = 8—1(x*,,, — x*). 


* 


. Set § = |x, 


. Follow Section 2 with the above p and §. Take §(x*;p+) =0. 





The points x, & x, + (s,At)v, are the “new quadrature points,” merely the result of applying the motion 
at each quadrature point. They will not necessarily be the quadature points of the new edge. In fact, the 
xX; need not be colinear if g > 2. As a result, the new edge is determined by best-fit, the solution to 


qd 
‘ ~ a\2 
am n 
agen: ARE YA) 


We can simplify our problem considerably with the following lemma: 


The least-squares best-fit plane of {xp} p-1 CR contains the centroid ¢ = + eet Xp. 





For a proof of this result, see the Appendix. We seek the “best one-dimensional subspace” spanned by the 
columns of A. To find this subspace, we call upon the singular value decomposition 


A= UV" 
where U € R?*? and unitary, V € R?*4 and unitary, and 5; € R?*? of the form 


O1 


xp 0, (q-r) 


y= =; = 


02-1) xr O(2-r) x (q—r) Z 
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in which r = rank A and 0, >--: >o, > 0. This second lemma further aids us in our quest: 


The first s columns of U represent the “biggest” s-dimensional subspace, the one that has the 
smallest representation error in Euclidean norm, contained in range(A). 





The proof of Lemma 2 is again contained in the Appendix. These two lemma explain the mysterious 
subtraction of ¢ from each column of A and c’ from A’. First, we note that range(A) = range[u; --- u,]: 


Ax = U|oiV?x --» opV?x O --- O 
= So ailx, vi)ui 
i=1 


Since the centroid is contained in the best-fit plane, by subtracting it from each data point creates a candidate 
for the best-fit line we seek. The task therefore becomes to find a single vector that best encapsulate the 
columnspaces of A (and similarly for A’), a problem to which we can directly apply Lemma 2. 


The above machinery explains how steps 1-5 of Algorithm 1 compute the new edges; we have € ++ (c,n;) 
and €’ ++ (c’,n2), where the pairing indicates the normal vector and point inside the new edge “plane.” 
From the contents of these pairings we can write parametric expressions for the new edges: 


Er+ Us) =c+s[ - n®), nr] 
E's U(t) =e’ +t[ - ne), ns] 
Step 6, therefore, computes the intersection of @ and @’ in terms of the above parameterization and sets the 


location of x*.y, to be the point ¢(s*) = é’(t*). As promised, the direction of motion p is the vector that 


new 
takes us from the original location to x{,,,; normalizing this vector and saving the original length gives us §. 


Since this completely describes the desired displacement, we quash movement in the orthogonal complement. 


motion 





Fig. 2: Illustration of Algorithm 1. The left shows an element with two receding sides (in red and 
blue) along with the quadrature points on each edge. The corner node x* is the purple star. Algorithm 1 
computes the receded edges independently and then finds the location of the new node by intersecting the 
lines formed by a best-fit plane through the perturbed quadrature points, depicted on the left. Note that 

recession need not be uniform across an edge, as shown here. 


While Figure 2 above depicts the case of a double receder, the algorithm cares not for the specific combination 
of sliding and receding edges. When one side is fixed, Algorithm 1 will not display the correct behavior when 
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the other side is receding; in implementation, it is advantageous to detect the presence of a fixed side set and 
immediately suppress all motion, eliminating the unnecessary overhead of the algorithm. When both sides 
are sliding, the end result is that the node will remain in its original location, so including logic to process 
that case quickly is also a boon to speed, although Algorithm 1 will also produce this effect. 


2. Interior Nodes: Receders 


This section will look very similar to the preceding one when written. 


Summary of method: Recede each edge by best-fit line as we did in the corner case. However, the intersection 
between these edges may not be defined, if they were parallel to begin with and remained so after recession. 
Instead, we use the internal edge from the element home to the receding edges. The intersection of the new 
edge with the internal edge is computed. This gives two new points in total; the centroid of the two is taken 
and that becomes the new node location. 


3. Interior Nodes: Sliders 


A slider is characterized by forbidden motion in the normal direction. Up to this point, we have avoided 
asking for this elusive vector at nodes. If the boundary were smooth, there would be no issue. As elements 
are polygonal, the boundary of the computational domain is only Lipschitz, which means that the normal is 
only defined almost everywhere. The set of points where it is not defined is precisely the set of nodes—and 
precisely where we need it. 


Because faces are straight lines between nodes, the normal is the same everywhere on a given face. On a 
particular face, the normal at the node can be defined via limiting procedure: for a fixed € € OOM OT: 


v(€é > x)= lim — v(x, + (6) — @1)h, v2 + (2 — 22)h’) (17) 


h—0+t ,h’>0- 


Unfortunately, it is insufficient to consider only a single element, as shown in Figure 3. 





Fig. 3: Difficulties in Finding v. In (a), for an “ordinary” 7 (shaded in green), (17) will always 
converge to the green vector independently of £1 € 07 N OQ. The steps of the limiting process are the 
dashed vectors. On J’ (shaded in purple), the purple vector will be obtained regardless of the starting &2. 
This results in a node possessing a normal vector for each element in which it is contained. 


On a per element basis, there is no issue in (a). However, the finite element basis functions a, are constructed 
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from shape functions; using the new basis functions {@,}, we’d have something like 


Per (x) xET 
Bele) =) Pure) xET! 


if the green node x* is node number k, within 7, kg within 7’, and k* overall. If different normals at x* 
are allowed for different elements, we would have ,, 7(x*) # z,,7/(x*) so that W,. would no longer be 
continuous. Consequently, the “ordinary case” of (a) must be resolved by simultaneously considering all 
elements containing each boundary node. Fortunately, in R?, there can be only two such elements. 


Let x € OT NOT’ and v1 = v(€z > x) and vp = v(E7, > x) be the normals obtained from (17) applied 
within 7 and 7’, respectively. We then define the normal at x to be the “average” of the two normals: 


V1 +V2 


v(x) (18) 





[vy + v9 


A visual depiction of (18) is shown below in Figure 4. 





Fig. 4: The Solution to v. v, (in green, solid) and v2 (in purple, solid) are produced from the limiting 
process (17); v (in red, dashed) is the result of (18) applied to these two vectors. 


D. Example Problems 


This section will contain several examples of the two-dimensional scheme in action. 


Ill. Three-Dimensional Meshes 


The straightforward analysis possible in the 2D case can be attributed to two key facts: 
1. Each node can be contained in at most two elements that intersect the boundary non-trivially. 


2. Given a direction of motion, there are exactly two vectors orthogonal to it. Choosing appropriately 
allows us to construct a symmetric transformation matrix M. 


Unfortunately, neither of these hold in three dimensions. Consequently, the script for analysis diverges 
significantly from that 2D. An example of how (1) above may be violated is depicted in the figure below. 


Another distinguishing feature of the 3D case is that there is more than one type of corner. In 2D, a node 
is a corner if it is either contained in two side sets or is contained in two faces within the same element. In 
3D, because the number of containing elements may be arbitrary, attempts to characterize nodes based on 
elements is an exercise in futility. Instead, nodes are classified according to the user-specified side sets: 
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(a) (b) 


Fig. 5: An Arbitrary Number of Boundary Faces: Consider the red node in each subfigure. 
Well-placed quadrilaterals lead to the “nicest” case of 4 containing elements, shown in (a). Addition of 
mixed elements, such as tetrahedrons and pyramids in (b), complicates things considerably. 





Definition (Node Classification). 





A boundary node is a pure corner node if it is contained in three side sets. It is an edge corner 
node if it is contained in two side sets. An interior node is contained in a single side set. 





(a) (b) (c) 


Fig. 6: Corners in 3D. Each side set is indicated by a different color. (a) depicts a pure corner. In this 
example, all three containing faces are in the same element, but this is not always the case. (b) and (c) 
depict edge corners. A bevy of configurations is possible, but there can be only two side sets involved. 


The use of side sets to classify nodes shifts burden onto the user during grid generation. This is necessary 
because there is no way to tell from the mesh geometry alone how a node should be treated. By using side 
sets for identification, control remains with the user so that no errenous processing will occur. 


A. General Framework 


As we did in 2D, we begin with the vector-valued shape functions on an element T 


wr (x) e1 ifl<k<N 
,(X) = Up-n(xX)eg if N+1<k<2N 
Ue—2N (X) €3 if 2N+1<k<3N 
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and define the new shape functions 





V(X) Pe l<k<m 
We(x) er m+t1l<k<N 
~ Wr—n(X) Gk—n N+1<k<m+N 
Pi (X) = 
Urn (X) €2 m+N+1<k<2N 
Wr—2N (X) Pe-2N I2N+1<k<m+2N 
Ur-on(x)eg  mt+2N+1<k<3N 





The directions pz, qx, and rz, should be considered arbitrary at the moment. We mirror the procedure that 
led to (6) and seek a matrix that encodes the transformation from old to new elements via 


3N 
Px(X) = S> Mauipy(x) 
l=1 


To extend the previous result (??) to three dimensions, we need to construct M to be symmetric. In 2D, 
restriction to a single direction of motion makes this task straightforward, as creating an orthogonal basis 
yielding a symmetric M is merely a matter of choosing the correct sign. In 3D, it is necessary to enforce 
multiple directions of motion. In such cases, a symmetric M is beyond the realm of possibility, as that would 
require that every invertible matrix be similar to a symmetric one. Consequently, developing the framework 
for obtaining the correction to the stiffness matrix Al, when M is not symmetric is our foremost task. 


B. Applying the Corrections 


1. Element Stiffness Matrix 


Suppose the directions of motion pz, qx, and rz, have been predetermined and that these directions are 
linearly independent. Furthermore, suppose that pv #0, gq? # 0, and 7?) # 0. Given a set of vectors that 
do not meet this criterion, we can always relabel them according to the following: 





Lemma 3 (Rearrangement Algorithm). 





Let Vk,, Vk., and vz, be linearly independent vectors in R?. Let 


(ki,k5,k4) € argmin  trace(diag Q,)7! (19) 
mE (ki ,k2,k3) 


where II(k1, k2,k3) denotes the set of permutations of {k1, k2,k3}, diag(-) is the operator that 
extracts diagonal elements from a matrix, and for a = (ki, kJ,kJ) the matrix Q, is 


ee? 


A 
= (1) (2) (3) 
” | Une | Ope VCR 


(3) 


2 
a leg 


(1) y?) 


| Uti | yp?) ’ 
% 
kg kg kg 


Then oy #0, Oe # 0, and Ue #0. 





15 of 30 


American Institute of Aeronautics and Astronautics 


The result can be seen clearly by writing out the trace (19) explicitly: 


1 1 1 
ly] 


trace(diag Q,,)7' = 
| ky Uns | Ug 


(20) 


If one of ol), Ors or oh is zero, (20) will be infinite, and thus cannot be the minimum unless 
A, 


all permutations result in this degeneracy. We now show that since span {vx,, Vk; Vk3} = R°, at 
least one arrangement results in finite (20). Let n,(v) count the number of zero components of a 
vector v € R?. WLOG, we may suppose that nz(vz,) > nz(VE.) > Nz(Veg)- 


Case I: nz(vx,) = 2. 
Because span {Vx,,Vk2, Vk3} = R®, we have 


det [vi Vio Vis | #0 


Let i be the nonzero component of v;, and 2’ amd 7” be the zero components. Then performing 
a Laplace expansion on the m* column to evaluate the above determinant, we have 


a 2 3 
Ub, det Ye ") | #0 
2 3 
=A 


or that det. A 4 0. Consequently, each column of A has at least one nonzero entry; hence, at least 
one of 7 = (ky, k2,k3) or w = (ki, kg, ko) makes (20) finite. 


Case II: nz(vz,) = 1. 
At least one of up # 0, ue # 0, or of # 0 holds; suppose one is mse #0. Set ki =m. It then 
follows that only one of Re) # 0 or Ro 0 is true. 

Subcase II.1: v?) = 0. 

At most one of ve) = =0Oor v2), = = 0 is true, where {m’, We {ki, ko, k3}\{m}. 


Subcase II.1.1: v?) # 0 and v?) # 0. At most one of vw) ) =0or re ) =0. Put 


ki! a m! i= az kl _ ml”! v&), = =0 
2 Lm" v®) =0 7 om v®) =0 


If y2)y@) 


Uy, #0, then we may assign m’ and m” to ky and ky a Sea 


Subcase II.1.2: stmt = = 0. Because nz(Vm') < 1, it must be that vw) # 0. Put 
kg =m’. Then ve? ) £0, so put kg =m" to give a permutation that makes (20) finite. 


Subcase II.2: v(3) = 0. Repeat Subcase II.1 with the symbol swap (2) (3). 
Case III: n,(vz,) = 0. 


All components of vz,, Vz,, and vz, are nonzero. All permutations yield a finite (20). 
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With the directions properly arranged by the Rearrangement Algorithm, we put 


pi Q|o -:- 0 po 0 0: 0 p> 0 0 -- 0 





(21) 











By the Rearrangement Lemma, each diagonal block of M is nonsingular. A secondary secondary benefit of 
assigning the directions in this manner is that the matrices 


py» va ry 
Mii = me M33 = me Ms; = 
nf) 2 rc 


should be reasonably well-conditioned: if one of Jo? <1, then per (20), we will have trace(diag Q,)~! > 1. 
Such a permutation is unlikely to yield the minumum trace, with one exception: a poor quality grid may 
yield a column in Q, with nothing but small entries. In this case, there will be no avoiding poor conditioning. 





Invertibility of M is necessary to return to the original finite element basis y, from ;,. First, observe that 
we may permute M into the following form: 





where Cy is the change of basis matrix 


PE 


A 
Cr =| qf 


T 
x 
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The directions pz, qx, and rz, are assumed to be linearly independent. Thus, each Cy, is invertible, and so 
is M’. As there exists a nonsingular matrix P such that M = PM’, it follows that M is also nonsingular. 


Partition M along the thick lines of (21) as follows: 





Mi,/0 Mi3 O Mis O 
0 I 0 0 0 0 

ji Ms; |0 M33 0 M35 0 alA B | 
0/0 0 I 0 0 ele | 
Ms:|0 Ms; 0 Ms; O 
0 0 0 0 0 I 


where the blocks of the rightmost side correspond to the dividing lines on the left side. By construction, 
A = M1, is invertible. To see that D~! exists, we permute M into the form 





I 0 0 0 0 0 
0 M33 O Msz35 O |] Ma, 

epee 0 O° Ee 10%, 0) 02. j) IP 
0 Ms; 0 Mss 0| Ms 5 
0 0 0 0 I 0 
0 Miz 0 Mis 0} My 


and confirm by direct computation (which we omit here) that M” has an LU-decomposition because M33 
and Ms; are invertible. We then recall that a matrix has an LU-decomposition if and only if all leading 
principle submatrices are nonsingular [Bernstein, p.260]. Because D is such for M”, it must be nonsingular. 


Recall the Schur complement S and “Haynsworth complement” H of M: 
S2A-BD'C 
H*=D-CA'B 
We can further compute the determinant of M by block: 
det(M) = det(A) det(H) 
= det(D) det(S) 
which gives that both S~! and H~+ exist. Knowing this, we recall the following result [B., p.44-45}: 





Lemma 4 (Block Inversion). 





Let AE Cr*", BeECcr™™, CEc™™”" andDEC™*™. If A and S are nonsingular, then 


so | -~S—1BD-! 








—D-'cS-! | D-!4+D-!cs-!BD-! 


If D and H are nonsingular, then 


-1 





A=! + ABH CA | —A*BH-? 





-H-'CA | a 
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Each My; € R™*™ is diagonal and I denotes the identity matrix of size N — m; the sizes of the zero blocks 
follow from context. To the return to stiffness matrix in terms of the old shape functions from the new, we 
must use equation (8) and not (??) because M is not symmetric in this case. As in 2D, we have 


Al; =MTKM 
=K+M 7PM! (22) 


Depending on the classification of the node, one or more directions may be unenforced. That is, we allow 
motion in that direction to be determined by the dynamics of the system. Because the order of directions is 
set by (19), P is dependent on its classification. Fortunately, for each boundary node x,, we can write the 
general form 


-1 T T T 
P,=e [foxx €,,€;, + fin ©;,%, + fa,k e:,€4,| 


where ig = k, i) 2K4+.N, and ip = k+2N, and 


Fe 1 px enforced ee 1 qx enforced Pi 1 ry enforced 
ae 0 otherwise ne 0 otherwise 0 otherwise 


so that the full form of P is given by 
P=. Pi (23) 
k=1 
Consequently, we need only compute 
= = / us T 
M~* (ej e7 M7! = (M~Te;) (M~*e;) 
for 7 = 10, 21,72 to generate everything we will need to assemble Al, on the fly. 


Ordinarily, inversion of a matrix is a tall order and recommended only if it can be done in exact arithmetic. 
Fortituous for us, because the blocks of M consist of entirely diagonal matrices, M77 is not nearly so 
daunting to construct. Observe that within the Block Inversion Lemma, the blocks are the same size so that 
S*=A7*+A "BH CA Sop SA BA: 
Ho = Des De CS= *BDo* DMS =H CA? 





as long as all these quantities exist. We previously established that they do; consequently, we may “pick 
and choose” which among these is most convenient to compute. Because S is m x m and diagonal, it is 
trivial to invert. Formation of H requires inverting A, again m x m, but then we must invert H, which is 
(N —m) x (N —m), the same size as D. Hence, there is no advantage offered by the alternate form of the 
Block Inversion Lemma. As D is already known, we choose to invert it. We begin by partitioning D as 




















I 0 0 0 0 
0] M33; 0 Msgs | 0 I; 0 |0 
D=|o] 0 I 0 /0}=] o0/D, 0 (24) 
0 Ms3 0 M:;5 0 0 0 I 
0 0 0 0 I 








so that we need only compute D;>!. This time, however, is is convenient to apply the “pick and choose” 
version of the Block Inversion Lemma. Doing so, we can obtain an expression for D>: 


s-! 0 —S~'M35M;,! 
D,* = 0 I 0 (25) 
—H-'M;3M3;, 0 H>! 
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where S is the Schur complement of D,. and H,. is the reduced Haynsworth complement of D,,., 


S = M33 — M35M;.'Ms3 
H,. = Ms5 — Ms3M33'M35 


With D~! in hand, we can write expressions for S and S~!BD=!: 
S = My, — Mi38~' (M31 — M35M5,'Ms1) — Mi5H,-! (Msi — Ms3M33 M31) 
spp =|0 -s-p, 0 -s-'P, | 
where 


P, =M,5H7'Ms3M33 — M,3S~* 
Py = M138~'M35M;. — Mi5H;! 


To get MT from here, it is simply a matter of multiplying everything out and transposing: 


Sa! - 6 S7'P3 ) SPs ) 

0 I 0 0 0 0 

(yi S'P, 0 S-1 4 S7!P3P35 0 S~'PyP;—H-'Ms3M3, 0 
0 0 0 I 0 0 

SP, 0 S~*P3P5-—S-'M35M;,' 0 H>-! + S-1PyP¢ ) 

0 860 0 0 0 I 


where 


P; = —S~* (M3; — M35MzMs1) 
Py = —H71 (Msi — Ms3M33 Ms1) 
Ps = Mi5H-'Ms3M3,' — Mi3S~! 
Po = Mi3S~'M35Mz. — M,5H7! 


(28) 


In implementation, due to the diagonality of all matrices invovled, full matrix multiplication may be avoided 
entirely by working node-by-node. To this end, we compute the ith diagonal entry of all the ingredients 


needed to create M7?: 

















git) = gg) — ere pei) — (1) _ Gere 
= 4; 3) 3 OG | % 7) 
(3) (2) (3) (1) 
He) —,8 4% pe) — _ 1 GQ) Mi 4 
neo a 40 Ga (" a) 
pid Pere py pen Pere Be 
( = 


= He?) SG He) 2) SG 


pan _ 24 9 pin _ 20 9 


2 SGi)\p ®t) 6 sGa78) AGO 


sage ewe ah Bae 
= DP; Sa) q; 7) He) i q?) 
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Noting that when a matrix multiplies e;, the result is the extraction of the k*" column, we have 





T 
Mole = [0...: 0, 8%, 05.40, B?),0,. 05.50, B® ,0,. 50 (33) 
ee 1 Nak 2 ios oye 
M-Te;, = eae (34) 
k-1 N- N-1 N-k T 
MTe;, = [0,.-.,0,5,9, 0,. "0,62 Dior g Ob Oe 0) (35) 
where k-1 N-k N-1 N-k 
k,k k,k 
GO a _ 2 a) _ Ps 
k Sh) Vk S(&k) k S(&k) 
32) = pe Cie he Berieye Sic, Bae ae ti. 
k ~ S(k,k) Vk ~~ O(k,k) S(kk) k ~~" (kk) HEH 2 
a) = Pa 0) _ PSP ge sa ere ee 
k ~ S(hR) Te ~~" Slik) Sth) -@) be S(hK) 


Finally, to assist in assembling (22), let 
Ry, = MTP, M! 


3 
= = = T 
=e 1S fia (M e;,) (M "ei,) 
j=l 
Per (34)—(35), there are only nine nonzero entries in Ry: 
RG) _ gol forBo Be a fie yey? ne fapon 5) 
Riots) = aot forB BC 4. Frwy ae fo,n0 6) 


Rio) = gl forbo Be a Fi, RYE (3) yi) ‘af fond 6) 


Riv) rz et fonBo pe ai fi nak (2) Ste fone 6) 
pve) _ et fon Be bop gy OAO spine? 62] (36) 


RO) - at fonB© Be 4. fi KV (3) y?) ae ford? 5) 














pias) ae fouBe pe bps. nyo (3) + fond 5) 
pis) ee fonB? Be A Fey? a foe 5) 


Ror) a gol fonB© p& mis fey ey i fo,n00 62) 


























Finally, recalling (23), the contribution to the global element stiffness matrix from 7 is then 


A|,=K+5—R, (37) 


2. Right-side Vector 


The final piece to the puzzle is the contribution to the right-side vector P| 3 for this, we need only compute 
the perturbation MB, where 


T 
B2| fox8(p1,..-;Pm)™ O fin8(Q1,---,dm)? O fox8(t1,--.,%m)? 0 (38) 
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Each 0 is a N — m row vector of zeros. The block 8(di,...,dm) € R™ is defined as 
§(x1; di) 
8(di,...,dm) + : 
3(Xm; dm) 
The inclusion of a direction d; indicates that the recession value §(x;) is in the direction d;. We denote it § 


instead of the usual s because this value may be nonphysical and encompasses all contributions to motion. 
For example, the recession of a quadrature point is sAt; § would include this factor of At. 


The multiplication i ) £ MB has the same block structure as (28) and (38); the result is 


T 








3Dc) 
Bon =| Bf o ES 0 Ef 0 (39) 
where the vectors E,, Ez, E3 have components 
EY? © fo 1.8(xi; Py)pO + fin8 (Xe: q;)q + fo,n8 (xis ry) 
EY = fo,n8(xi3 Py)p? + F165 (5 q;)q\” Tr fo,n8(xi3 ry)? (40) 
By) = fo,n8(xi3 Py pO + fi,n3 (xa; q;)q a fon (xi3 ry)r 


3. Summary 


For node & within an element, suppose we have determined the directions of motion v,%,, Vz,, and vz, with 
known recession values §;, §2, and $3, respectively. These values and directions may be nonphysical and may 
not all be enforced. The procedure to compute the needed corrections is as follows: 

1. Arrange the directions into pz, qx, and rz via the Rearrangement Lemma. 

2. Mark which of these directions are to be enforced via fox, fi,x, and fox. 


3. For each node on the boundary xz, 1 < k <_m, compute the nonzero entries to the correction matrix 
Rx via (36). Add these values into the local element stiffness matrix A|_- 


4. Compute the correction (39) to the right-side vector. The entries of this addition are given in (40). 


C. Quantifying the Motion 


Section B assumes that all the directions pz, qx, and rz and the corresponding § values are known a priori. 
CHAR specifies the values of 5 only at the face quadrature points. Constructing the direction and recession 
pairs (px, (xx; Px)), (de, (Xx; ax)), and (re, 8(xz;7%)) at each node x; from § and v at quadrature points 
is non-trivial and is the focus of this section. To illustrate the challenges we face, consider Figure 7 below. 





Fig. 7: Physically Correct Recession. If the temperature over the parallelogram’s surface (red) is 
uniform, the surface should recede in a direction parallel to its current plane of existence (right). 
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Figure 7 depicts a parallelepiped with a single receding face (colored red); the faces adjacent to it are 
sliding faces and the remaining one (parallelogram in the back) is fixed. On physical grounds, we expect the 
parallelogram face to melt away and the element to kept its general shape through recession. 


Figure 8 above highlights the difficulty in choosing directions of recession that lead to the physically correct 
result. At the edge corners (in blue), we require that the node remain on the sliding face; this keeps the 
exoskeleton formed by the four corner nodes and preserves the parallelepiped’s general shape. 


recede 


area 





Fig. 8: When the “Right” Direction Goes Wrong. Edge corners (blue nodes) should move in the 
face in the sliding side set so that the sliding condition trumps the receding one. At the red interior nodes, 
the “natural” direction of motion is the normal vector. The result of this motion is on the left. 


At the interior nodes (in red), we have no information of the existence of a sliding face because all analysis 
is done on a per-element basis. Furthermore, even if we did know of the sliding side sets, there is one in each 
cardinal direction—which is the “right one”? As such, it seems that the natural direction for interior nodes 
is the normal vector. However, the right side of the figure reveals that such motion causes mesh quality to 
deteriorate. Continued motion along the normal sets up a future unavoidable collision with an edge. 


Complicating matters is that each interior node is contained in multiple elements. Figure 8 depicts a simple 
example in which each face in the receding side set has the same normal. If the receding set is a surface and 
not a mere plane, defining the normal at nodes is difficult—and defining § is harder yet. 


1. Pure Corner Nodes 


Perhaps contrary to expectations, the simplest case is the pure corner. A typical situation is depicted below: 





Fig. 9: Typical Pure Corner Case. The three different colors indicate the three side sets that contain 
the pure corner x*. The normal vectors for each side are indicated. 


If any of the three side sets is of motion type “stationary,” the node should remain fixed in place. We set 
Y y yy 


P = €1, q = eg, r = eg and §(x*; -) = 0 and move on. We therefore focus on combinations of sliding and 
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receding side sets. Unlike the other node types we shall explore later, there are no subcases for each possible 
combination of sliding and receding side sets. 


It is tempting to choose the face normals as directions of motion (as in Figure 9). However, § is known only 
at quadrature points; choosing the value of s at the closest quadrature point is physically incorrect and likely 
to lead to poor mesh quality and skewed geometry. Instead, the solution is, in essence, to move entire side 
sets: we collect the face quadrature points from each face containing x* in the side set, move each according 
to its known values of v and s, and then compute the best-fit plane through the resultant points. We then 
find the intersection of these three planes, guaranteed to be unique except in degenerate cases. 


The process is prescribed precisely in Algorithm 2 below: 





Algorithm 2 (Pure Corner Motion). 





. Choose a side set S; which contains x*. Let 1,...,7m, € Sj; be the faces containing x*. 


. Let Xp1,---,Xk,q, € Fx be the quadrature points of the k*® face. Let vy; and &;; be the 
normal and recession value, respectively, specified at the quadrature point x, j. 
mj; 
. Compute the centroid of the collection U {Xi,i + (Sn At)V Ki 
k=1 


qk 
i=1- 


qk 
[xn,e + ($n,cAt)vr,i] 


. For each k = 1,...,m,;, form the 3 x gq, matrix 


Ax 7 XK1 + (Sp 1 At)vE — Cj | tee | Xk,qy + (Sk,qx At)V kay — Cj 


and glue all of these together to form the 3 x Q; matrix, where Q; 4 Pe dk 


A= A, «+: An, 
. Compute the singular value decomposition of A: 
A=UdV" 
where U € R®*3, © € R°*@i, and V € R@i*2s. Let d; £ uz, the 3™ column of U. 


. Repeat Steps 1-5 for the remaining two side sets to obtain cz, de, cs, and d3. 


. Solve the linear system 
df (di, ci) 
d3 Xnew = (da, C2) 


dz (d3, ¢3) 


=| 


. Set A= x*,.. — x* and § 4 |A|. Compute the singular value decomposition 


new 
a=UvV" 
A ads “ Aa Aw A oid th a 
. Set wi = (sgn(a, f1)) G1, wo = Ge, ws = tis, where i; is the i‘ column of U. 


. Follow Section B with the direction and recession pairs (w1, 8), (w2,0), (w3,0), all enforced. 
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Steps 1 and 2 are self-explanatory: collect all faces containing x* along with the quadrature points and 
recession values and normals at those points. In steps 3 and 4, we work with the “new quadrature” points 
Xk 4 Xr + ($x iAt)vn;. Because $ is specified at each quadature point, we can compute with absolute 
certainty the new locations of only quadrature points after recession. 


From the collection of X,,; from the entire side set (containing x*), we must extract a single plane. For all 
but sliding side sets or single face receding sets and those with uniform temperature (so that s = constant 
across the side set), the set of X;,; are unlikely to be coplanar. To determine the post-recession face, we 
must obtain a best-fit plane through these points. As in 2D, Lemma 1 and Lemma 2 aid us in our question. 
As we did before, we call upon the singular value decomposition to obtain this best-fit plane. 


Each plane generated in Step 5 is described by the equation (x —, d;) = 0. Writing out these three dot 
products and casting into matrix form leads to the 41, the solvability of which is aided by the following: 


The normal vectors v1, V2, and v3 of adjacent faces of a tetrahedron span R?. 





Let a, b, and c be the edge vectors at the vertex formed by the intersection of the faces having 
normals v1, v2, and v3, arranged such that vy; =a x b, v2 = bx c, and v3 =ax c. Because all 
vectors lie in R?, we have that 


det [ vy | v2 | vg ] =v1- (v2 x vs) (42) 


Expanding v2 and v3 in terms of a, b, and c, and applying the identity for the the vector triple 
product a’ x (b’ x c’) = b’(a’- c’) — c’(a’- b’), we have 


V2 X v3 = (a-(bxc))c 


Hence, 


V1 + (V2 X v3) = [(ax b)-¢] |(b xc) - al = [(ax b)-c]” 


by the nature of the scalar triple product. The volume of the tetrahedron is 2 |(a x b) - cl, so if 


the lefthand side is zero, then the tetrahedron is vacuous—a contradiction. Accordingly, it follows 
that v1, v2, and v3 are linearly independent and hence span R°. a 





For hexahedra and pyramids, we can construct a tetrahedron contained in the original polygon that overlaps 
with the original adjacent faces in question, thus extending the above result to these polygons. The lemma 
itself, of course, is not directly applicable to (41). The vectors that determine (41) are normal vectors from 
planes not known a priori to form a polygon. However, since the new plane is a best-fit of perturbed points 
from the original and the determinant is a continuous map that measures volume, we can reasonably expect 
(42) to remain nonzero. The contrary requires motion of one face so dramatic that it becomes parallel to 
another, highly unlikely unless this were nearly true originally (e.g., a shield-like corner). 


The intersection point determined in Step 7 is the location of the new node. Step 8 computes an artificial 
direction of motion a, directed toward the new location from the old. The recession value § is the length of 


: 3 : : ts ; ayt : : 
this vector; SVD provides both a normalized version of A and a basis for {a} . We must adjust for sign, 





since fy = + |aj~* a. It is crucial that the direction be pointed toward Xnew to be consistent with §. 


2. Edge Corner Nodes 


Recall that an edge corner node is one that is contained in two different side sets. Geometry can differ wildly 
among instances of edge corners (two simple examples are depicted in Figure 6). Since the set of potential 
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motion 


== 





Fig. 10: Visual Depiction of Algorithm 2. This figure depicts a simple example of Algorithm 2 in 
action. The geometry before motion is shown on the left. The blue (left) and yellow (right) faces are 
receding, with uniform recession at each face’s quadrature points but different between the two faces. The 
green (top) face is sliding. The left diagram shows the motion of the node: it moves from the purple circle 
to the red star Xnew. The direction A is indicated. The length of this vector gives §. 


configurations is seemingly unbounded, any attempt to account for special geometries amounts to a feral 
anatidae pursuit (or more colloquially, a “wild goose chase”). Instead, we design for “typical” cases such as 
those in the figure and place some burden on the user to have designed a reasonable mesh. 


We identify three possible configurations: when both side sets are receding, when one is sliding and the other 
receding, and when both are sliding. When one set is stationary, we set all motion to zero as we did before. 


3. Edge Corner: Receder/Receder 


We begin with the most complex of the cases: double receders. This configuration causes the most issues for 
simple mesh motion schemes; success in handling this case is therefore a hallmark of the tangle-free mesh 
motion scheme. The method is ostensibly similar to pure corner case but necessarily differs in a few aspects 
due to the increased complexity of the edge corner. 


Summary of method: Very similar in spirit to the pure corner case. We have two sides in motion. We 
compute the new faces by best-fit planes as we did in the pure corner case. To replace the missing third 
side, we use internal (non-boundary) faces that contain the node. The intersection with the new faces is 
computed for each internal face. We then take the new location of the node to be the centroid of all these 
intersection points. 


4. Edge Corner Nodes: Receder/Slider 


We begin by presenting the mathematical algorithm. Explanation of the various steps follows after. 





Algorithm 3 (Edge Corner: Double Receder). 





1. Let F1,...,Fm be the boundary faces containing x* in the receding side set; F{,...,F/,,, 
in the sliding side set. 


2. (Receding side set.) Let Xp1,..-,Xk,q, € Fx be the quadrature points of the k*® receding 
face. Let vy; and 8; be the normal and recession value, respectively, specified at x,.j. 
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continued). 


3. Compute the centroid of the collection U7", {Xx,i + (aecAt)vna he 
1 m dk 
ct iss S° yD [Xk ,i + (S:At)vi.,| 
k=1 w=1 
4. For each k = 1,...,m, form the 3 x q, matrix 
Ax 4 Xk + ($, 1 At)VvE 1 —c | tee | Xk,q, + (Sk.q, At) Vk, —c 
and glue all of these together to form the 3 x Q matrix, where Q = eer dk, 


A=[A gts An | 


. Compute the singular value decomposition of A: 
A=UdV" 


where U € R°*3, © € R°*®@, and V € R@*®. Let P be the plane with normal vector 
v = us, the 3° column of U, passing through the point c. 


. (Sliding side set.) Let x), 1,.-. Kh gf € Fi, be the quadrature points of the k‘ sliding face. 


. Compute the centroid of the collection Up_, {x4} 75: 


. Form the 3 x Q’ matrix, where Q! = ga Gus 


/ / 


rh / / 1 
A’ = X11 —c hei ty | X1 gy —c wei | Xai! ra 
m 


. Compute the singular value decomposition of A’: 
A=U'S'v"™ 


where U’ € R°%3, SY € R3*@", and V’ € R@’*@’, Let P’ be the plane with normal vector 
vy’ =u, the 3°? column of U’, passing through the point c’. 


. Let Fy/,...,F;/ be the non-boundary (internal) faces containing x*. Let vf and ci? be the 
normal and centroid, respectively, of Fy’. 


. For each k = 1,...,b, compute the intersection of P, P’, and F;! by solving 
7 [me 
oe esle.a 


(UK Cx) 


new’ 


b 
. Compute the centroid of {xi2 . Call this point x* 
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continued). 


13. Set A= x* 


* ww — X* and § = |A|. Compute the singular value decomposition 


a=Usv" 


th column of U. 


A A A A AA Aa A : . 
14. Set wi = (sgn(a, i1)) G1, we = Go, w3 = tis, where ii; is the 7 


15. Follow Section B with the direction and recession pairs (w1, 8), (w2,0), (w3,0), all enforced. 





Steps 1-5 mirror those for the pure corner case exactly; steps 6-9 repeat the procedure for the sliding side 
set. As before, we compute a single best fit plane for all containing faces in each side set. In the pure corner 
case, Here, though, the maneuver seems to be more precarious. 

[Some results and theorems that prove what we did isn’t so bad] 


5. Edge Corner Nodes: Slider/Slider 


This section will look similar to the one on Pure Corners once it is written. 


Summary of method: The two edges, €, and €2, containing the node that live between the sidesets are 
identified (there are always only two). We move along €; to the next node, identify the other edge between 
the sidesets containing that node, and move along it. We repeat this process until we come to another sideset 
or loop back to where we started. If the encountered sideset is a receder, then we select €2 as the direction 
of motion; otherwise, we use €;. We then compute the orthogonal complement to the selected direction with 
the SVD and enforce zero-motion conditions in these two orthogonal directions. 


6. Interior Nodes 
Recall that an interior node is defined by a single common side set shared among all faces containing x”. 


Our analysis requires us to distinguish three cases by the nature of motion of this side set. 


The easiest case is when the motion of the side set is “stationary.” While this case can be covered by one 
of the two coming cases, the amount of detection and computation done in those cases ultimately goes to 
waste for a stationary side set. It is therefore more efficient to detect this case and handle it immediately: 
set p = €1, q = €2, r = e3 and §(x*; -) = 0 and enforce these boundary conditions. 


7. Interior Nodes: Receders 


It is most clear to present the algorithm first and then explain its components afterwards. 





Algorithm 4 (Interior Receder). 





. Let Fy,,...,Fm be the boundary faces containing x*. 


. Let Xp1,---,Xk,q, € Fx be the quadrature points of the k*® face. Let vy; and &;; be the 
normal and recession value, respectively, specified at the quadrature point x; ;. 


dk 


. Compute the centroid of the points {xxi + (8, ,At)V xi coal 


1 qk 
Cr es ac ‘> [ki - ($2,.At)vx,,| 
* i=1 
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continued). 


4. Form the 3 x q, matrix 


Ay= Xei + (Spc At), — Ce | mee | Xkax + ($k,q, 46) kg, — Ck 


. Compute the singular value decomposition of Ax: 
A; = U,54 Vi 
where U;, € R®*3, 5, € R?%%, and V;, € R®*%. Let &,, * uz, the 3° column of Ux. 
. Repeat Steps 2-5 for each k € {1, hi ,m} to obtain c;, and vp; for each face Fy. 


. Locate each edge containing x* that is not contained on F;,. This will yield the collection 
{ex i}, each of which should be first normalized. 


. For each i = 1,...,,, compute the intersection of the line ¢;(t) = x*+t€,; with the plane 
with normal %; passing through the point cz (k* best fit plane). First, compute 


A (cx — x*, D;) 


a, = = 
‘ (Ex,i, ae) 
The intersection point is then given by 


A * = 
Yri=X + Ue 


. Repeat Steps 7 and 8 for each k € {1,...,mb}. 


m 
. Compute the centroid of the collection U 1 Vice be 
k=1 


a —-l1 m “mn, 

x oA 

Xnew — Nk Yk,i 
k=1 


k=1 i=1 


4 |A|. Compute the singular value decomposition 


a A Aa 
. Set A= x%.,, — x* and § 


a-UsvVT 


th 


A A A A Aa Aa A : . . 
. Set w, = (sgn(a, f)) G,, w2 = Ge, w3 = tz, where fi; is the i" column of U. 


. Follow Section B with the direction and recession pairs (w1, ), (w2,0), (w3,0). 





8. Interior Nodes: Sliders 


This section will look similar to the preceding one once it’s written. 


Summary of the method: use radial basis functions to obtain an estimate of the surface normal vector. By 
perturbing along the known normal vectors at the quadrature points, we can construct a function whose 
gradient gives the normal at any point. 


D. Example Problems 


This section will contain several examples of the three-dimensional scheme in action. 
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Fig. 11: Algorithm 4 in Action, Part I. This figure shows Steps 1-6 of Algorithm 4 in action. The 
quadrature points on each face move their specified s units along the normal for that location. The best-fit 
plane for the new points is then computed (colored planes) by the SVD. 


IV. Conclusion 


Summary and wrap-up. 


Appendix 


Some proofs of results I have used; they are not new results but the insight in the proof may be useful to 
those who have not seen these concepts before but are interested. 
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